Local Density Functional Theory for Superfluid Fermionic Systems: The Unitary Gas 
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The first detailed comparison between ab initio calculations of finite fermionic superfluid systems, 
performed recently by Chang and Bertsch [Phys. Rev. A 76, 02f603(R), (2007)] and by von Stecher, 
Greene and Blume [e- print arXiv:00705:067fvf], and the extension of the density functional theory 
Superfluid Local Density Approximation (SLDA) is presented. It is shown that SLDA reproduces 
the total energies, number density distributions in inhomogeneous systems along with the energy 
of the normal state in homogeneous systems. Unlike the Kohn-Sham LDA, in SLDA the effective 
fermion mass differs from the bare fermion mass and the spectrum of elementary excitations is also 
reproduced. 

PACS numbers: 31.15.Ew, 71.15.Mb, 03.75.Ss 
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The Density Functional Theory (DFT) introduced by 
Hohenberg and Kohn [l[ became the tool of choice in the 
calculation of the properties of essentially most electron 
systems Q after the introduction of the Local Density 
Approximation (LDA) by Kohn and Sham In order 
to achieve the accuracy needed in particular in chem- 
ical applications, a number of extensions of the LDA 
have been developed, the Local Spin Density Approxi- 
mation (LSD (A)), the Generalized Gradient Approxima- 
tion (GGA), etc., which have been thoroughly tested on 
a large variety of systems over the years by comparing 
the results of the LDA, LSD(A), and GGA with ab ini- 
tio calculations and by refining the form of the density 
functional used in practice 0, 0| • All of these formu- 
lations rely on the Kohn-Sham orbitals and thus cannot 
deal effectively with superfluidity. The DFT extension to 
superfluid systems is a fundamental problem of the many- 
body theory. Almost two decades ago such an extension 
was suggested Q , however in terms of a non-local pairing 
field. A DFT formalism in terms of non-local fields is def- 
initely intuitively less transparent, significantly harder to 
deal with in practice, and most likely physically not well 
motivated. The fact that the BCS theory leads formally 
to a non-local pairing field is not a strong argument in fa- 
vor of such an approach. Such an argument was not used 
in the case of normal systems. All the BCS results could 
be recovered easily within a fully local formulation 0, 0] • 
There was a technical motivation to proceed with such a 
non-local formulation in the case of superfluid systems: 
the presence of a rather annoying divergence in the defini- 
tion of a local anomalous density [1, 0] , an issue which has 
been successfully dealt with in Refs. @, 0, HI • This par- 
ticular extension of the DFT to superfluid systems was 
dubbed Superfluid Local Density Approximation (SLDA) 
and it was applied so far to nuclear binding systematics 
and to quantized vortices Q . 

Here I shall analyze the properties of a fermion system 
at unitarity, when the strength of the interaction in a 
two-species fermion system (spin-up and spin-down) cor- 
responds to an infinite scattering length [9(. The basic 



properties of such a system, the energy per particle, pair- 
ing gap, etc. have been established in a series of ab initio 
calculations^ EH, El El for homogeneous systems. In 



Refs. [14j, [15[ two independent groups report on ab ini- 



tio calculations of the properties of fermions at unitarity 
in a harmonic oscillator trap. These results provide the 
unique opportunity to test directly the accuracy of the 
SLDA, by deducing from the ab initio calculations of the 
homogeneous matter the appropriate local energy den- 
sity functional and then use it to predict the properties 
of a fermion system in the presence of an external one- 
body potential. The proof of the Hohenberg and Kohn 
theorem is extended in a trivial manner 2] to the case 
when the external field has in second quantization the 
structure: 

]T K, t (r)^(r)^(r) + [A ext (v)yj\(v)yj\(r) + H.c.}. 

(1) 

Thus one establishes that there is a unique mapping be- 
tween the external potential, the total wave function of 
the system and the normal and anomalous densities and 
that a unique density functional of these densities ex- 
ists. Consequently, the ground state energy of a super- 
fluid system can be computed using a functional of the 
normal and anomalous densities. Here I shall analyze 
systems for which only V ex t(v) exists and one can con- 
sider formally that the external pairing field A ext (r) is 
vanishingly small, but not identically zero. This is as- 
sumed simply to force the reader to appreciate the fact 
that particle projection to a "wave function with a well 
defined particle number" is neither required nor needed 
in order to recover the correct ground state energy. 

Fermions in the unitary regime (when the scattering 
length a is large) are particularly attractive for a num- 
ber of reasons: (1) this is a strongly interacting system 
which exhibits superfluid behavior and a complex phase 
diagram [HI ]; (2) at unitarity (\a\ = oo) the form of the 
energy density functional is restricted by dimensional ar- 
guments; (3) the availability of ab initio results for homo- 
geneous and inhomogeneous systems; (4) the relevance of 
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this systems to a large variety of physical systems (low 
density neutron matter in neutron stars, cold fermionic 
atoms in traps, high-T c superconductivity, etc.); (5) tun- 
ability of the interaction both theoretically and experi- 
mentally. 

Dimensional arguments suggest at unitarity the sim- 
plest SLDA energy density functional (units h = m = 1): 
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where a, /3 and 7 are dimensionless parameters and 
n(r), r(r) and i^(r) are the number /normal, kinetic and 
anomalous densities expressed through the usual Bo- 
goliubov quasi-particle wave functions [wfc(r), Vk(j)] and 
where k labels the quasi-particle states. The new uni- 
versal parameter a is being introduced in this work and 
its presence proves to be critical in order to achieve the 
high accuracy demonstrated below. Since the kinetic and 
anomalous densities diverge @, 0] one has to introduce 
a renormalization procedure for the pairing gap and for 
the energy density. The renormalized density functional 
and the equations for the quasi-particle wave functions 
obtained by the standard variation are as follows: 
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+ A c (r), 
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(5) 
(6) 

«:(r)« fc (r), (7) 
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[h(r) - fi]u k (r) + A(r)v k (r) = E k u k (r), 
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and thus the pairing part of the functional is uniquely 
defined. The second term in Eq. (fT0"| for U(r) = 
5£(r)/6n(r) is obtained by varying g e ff(r), see Eq. ((6|), 
with respect to n(r) and neglecting in the first approxi- 
mation the dependence of A c (r) on U(r). There is still a 
small correction term to U (r), arising from varying A c (r), 
which can be made rather small if E c is sufficiently large. 
This additional term can be computed using 
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(15) 



and one can argue that to a good approximation 
5U(r)/5n(r) w 2(C/(r) - F ext (r))/3ri(r), which is con- 
sistent with the universality of a homogeneous unitary 
Fermi gas (UFG), where U(r) cx n 2 ^ 3 (r). The results are 
independent of the value of the cutoff E c if this is cho- 
sen appropriately large [l8j]. The above formulas apply 
to systems with an even particle number. In order to 
describe systems with an odd particle number one has to 
place an extra quasiparticle in a specific quantum state 
n [3. 




+ U(r). (13) 



FIG. 1: (Color online) The ab initio fermionic quasiparticle 
spectrum determined in Ref. [ill ] ( red circles) compared with 
the SLDA spectrum defined by the ab initio parameters £, rj 
from Ref. GJ] (solid blue line), [l(J (dotted blue line), £2 
(dashed blue line) and naive BCS approximation at unitarity 
(black dotted-dashed line). 



Here E c is an energy cutoff and I have added as well an 
external potential V ex t(j). It can be shown that in the 
total energy the kinetic and anomalous densities have to 
enter in the combination 



a- 



To(T) 



A(r)v c (r) = a 



r c (r) 



■ 5e# (r> c (r)| 2 (14) 



By requiring that a homogeneous gas of number den- 
sity n — N/V = kp/3Tr 2 has an energy per particle 
E/N = 3£s£_f/5, a chemical potential fi — £,s£f and 
a pairing gap A = r/eF, where ef = kp/2, one can de- 
termine the dimensionless parameters a, (3 and 7 in Eq. 
([5l). The corresponding equations determining these pa- 
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rameters are: 



kp 



d 3 k 
J2nf 



r d 3 k 




J (2tt)3 


I s ? ( 




r d 3 fc 




(2tt) 3 



where 



E k 



1 

Qffc 2 



A 2 ' 
2E k 

2E~J ' 



(16) 



-ej.n/3, (17) 
5 



(18) 



ak 2 
- (d 



U-fi-- 

A2^/2 



ak 2 



+ (0-Zs)e F) (19) 



(3tt 2 ) 2 / 3 77 2 
67 



(20) 



The momentum hp can be factored out of the equa- 
tions [TBI [17]), and ([18] . as expected for a system at 
unitarity. Using £ s = 0.42(2) and 7/ = 0.504(24) deter- 
mined by Carlson and Reddy [ll[ one obtains a = 1.14, 
(3 = —0.553 and 1/7 = —0.0906. Using the original 
values of Carlson et ai, [l(|, namely £s = 0.44 and 
i] = 0.486 one obtains instead a = 1.12, ft = —0.520 
and I/7 = —0.0955. A different set of values has been 
calculated recently by Juillet 0, £<? = 0.449(9) and 
7] = 0.442(3), which lead to a = 0.812, (3 = -0.172 
and I/7 = —0.0705. The value of a = 1.14 leads to an 
effective mass of vn e g jm = 1/a = 0.877. 

One can now calculate the spectrum of the elementary 
fermionic excitations of a homogeneous UFG, shown in 
Fig. [TJ and compare it with the ab initio spectrum deter- 
mined in Ref. ll| . The agreement between the two spec- 
tra is nothing short of spectacular, especially if one has in 
mind that DFT is not usually expected to reproduce the 
single- particle spectrum 0, 0, 0| , unless one sets up from 
the outset to achieve that as well 01. The value of 
the effective mass is a matter of mild surprise, and also 
the fact that the minimum of E k occurs at fco = 0.906kp. 
The analysis of Ref. 2l| arrives at a similar value for 
fco. The energy per particle of a homogeneous UFG 
in the normal state has also been determined 0, Il3j |. 
E N /N = 3£nSf/5 with £ N — 0.55, which is an average 
of the two results. The SLDA energy density functional 
is consistent with this value, as £jv = ol + (3 — 0.59. The 
parameters £ and from Ref. [l2j lead to £y = 0.64 
and an effective mass m e ff/m = 1/a = 1.23 and to a 
quasi-particle spectrum in strong disagreement with the 
results of Ref. [L^. 

The SLDA calculations for the finite systems [with 
V ex t(v) — muj 2 r 2 /2] have been performed using the 
Bessel Discrete Variable Representation (DVR) method 
fl7j |. and a few details are presented in Ref. [1S\. The 
comparison between the Green Function Monte Carlo 
(GFMC)_Jlil, Fixed Node-Diffucion Monte Carlo (FN- 
DMC) [ljand SLDA is presented in Figs. HH and 
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FIG. 2: (Color online) The comparison between the GFMC 
[3, FN-DMC [3 and SLDA total energies E(N). The 
clear odd-even staggering of the energies is due to the on- 
set of the pairing correlations. The inset shows the discrep- 
ancy between the GFMC and FN-DMC and SLDA energies, 
SE(N) = Emc(N)/E S lda(N) - 1, where E M c(N) stands 
for the energies obtained in GFMC or FN-DMC respectively. 



18| (units for energy and length defined by 
u> = 1). The agreement between the Monte- 



in Ref. 

h = m 

Carlo (MC) and SLDA results is very good, especially 
keeping in mind that the MC calculations for both infi- 
nite matter and finite systems have aside from statistical 
errors also noticeable systematic errors. Both GFMC 
14| and FN-DMC 15j calculations are in principle vari- 
ational and, since the energies for the larger systems in 
the FN-DMC calculation are consistently lower that the 
corresponding GFMC results, one can expect that the 
FN-DMC results are somewhat more accurate. That is 
also in line with the smooth behavior of the 5E(N) = 
Efn-dmc(N)/ Eslda{N) — 1 with N as opposed to the 
behavior of SE(N) = E GFM c(N)/E SLDA (N) - 1. The 
SLDA energies converge unexpectedly fast towards the 
FN-DMC values, see Fig. Hand Table I pj|- As a refer- 
ence, the Thomas- Fermi energy E(N) = (3A r ) 4/3 v / S : /4, 
which for N = 20 comes to 38.05 as compared to 43.2, 
41.35 and 41.51 in GFMC, FN-DMC, and SLDA respec- 
tively. This discrepancy, whose size is typical in this par- 
ticle range, is comparable in magnitude with the conden- 
sation energy, but opposite in sign. 

As far as the density profiles are concerned, they agree 
reasonably well in the surface re gion , but show noticeable 
differences in the central region [18| . The reasons for this 
disagreement are likely the relatively poorer quality of 
the GFMC results for larger particle numbers [lij]. It 
is notable that there is a small discrepancy between the 
SLDA energies calculated as expectation values of the 
functional and using the viral theorem ^22], a discrep- 
ancy which seems to decrease with increasing E c . No- 
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tice that there is a somewhat bigger discrepancy in the 
GFMC values of the energies and the corresponding virial 
expectations 14, Til]. 




FIG. 3: (Color online) The pairing field A(r) for even particle 
number systems. The inset shows the quantity 52E(N) = 
E(N) - (E(N + 1) + E(N - l))/2 calculated in SLDA (red 
squares) and GFMC (blue circles). 

Even though the agreement between the SLDA and 
ab initio results is surprisingly good, a better agreement 
would have been bad news. The reasons for the dis- 
crepancies can be ascribed to several origins: (i) The 
energy density functional assumed here, see Eq. ([5]), is 
certainly not unique and further efforts should be de- 
voted to study other possible forms, (ii) The so called 
Self-Interaction Correction (SIC) is not present here, 
and its absence is seen in the SLDA energy for N = 1, 
which is 1.37fku instead of 1.5huj. Even though S2E(N) = 
E(N) - [E(N + 1) + E{N - l)]/2 calculated in the two 
methods agree within the error bars, some differences are 
likely due to (Hi) the absence of the polarization of the 
even core by the field of extra odd particle (full spheri- 
cal symmetry was assumed). Another reason is (iv) the 
absence of spin number densities in this formulation of 
the SLDA and an extension to Superfluid LSD approxi- 
mation is needed, especially for small systems. For now, 
the energies of the homogeneous asymmetric systems are 
known with significantly less accuracy, see Ref. [lij ]. 
And last, but not least, (v) gradient correction terms 
are needed (beyond those already included through the 
use of explicit single-particle kinetic energy) . 

The differences between the ab initio and the SLDA 
results are noticeably less than the corresponding dif- 
ferences in electronic systems and the validation of the 
DFT extension presented here is to my knowledge a first 
of such kind. I anticipate that the existence of an accu- 
rate SLDA will have a great impact on several physics 
subfields: atomic nuclei, neutron matter and quark mat- 



ter in neutron stars, dilute atomic Fermi gases, condensed 
matter and, maybe, on high T c superconductivity as well. 
The ability to perform high accuracy calculations using 
essentially meanfield techniques in lieu of ab initio calcu- 
lations can hardly be underestimated. 

I am grateful to S.Y. Chang and G.F. Bertsch for the 
use of their unpublished results and for discussions, J. 
Carlson for the data from Ref. llj, to M.M. Forbes for 
numerous discussions and help with coding and numerics. 
A number of discussions and the valuable comments of 
W. Kohn are particularly appreciated and I thank D. 
Blume for information about their unpublished results. 
Support is acknowledged from the DOE under grants DE- 
FG02-97ER41014 and DE-FC02-07ER41457. 

Notes added in proof. The recent FN-DMC calcu- 
lation of D. Blume, J. von Stecher, and C.H. Greene, 
arXiv:0708:2734vl, extended to both even and odd sys- 
tems up to N — 30, demonstrated very nice agree- 
ment with SLDA. The Galilean invariance of Eqs. (2) 
and (5) becomes manifest upon adding the term (1 — 
a)p 2 (r)/2n(r), where p(r) = Real[— i^2 k u£(r)Vufc(r)]. 



[9 



[10 
[11 

[12 
[13 

[14 

[15 
[16 

[17] 



P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 
(1964); W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). 
R.M. Dreizler and E.K.U. Gross, Density Functional 
Theory: An Approach to the Quantum Many-Body Prob- 
lem, (Springer, Berlin, 1990). 

W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965). 

J.P Perdew et al, Phys. Rev. Lett. 77, 3865 (1996); J.P 

Perdew et al, Phys. Rev. Lett. 82, 2544 (1999); M. Seidl 

et al, Phys. Rev. Lett. 84, 5070 (2000). 

L.N. Oliveira et al, Phys. Rev. Lett. 60, 2430 (1988); S. 

Kurth et al, Phys. Rev. Lett. 83, 2628 (1999); E.K.U. 

Gross et al, AIP Conf. Proc. 577, 177 (2001). 

A. Bulgac and Y. Yu, Int. J. Mod. Phys. E 13, 147 

(2004) . 

A. Bulgac and Y. Yu, Phys. Rev. Lett. 88, 042504 (2002); 
A. Bulgac, Phys. Rev. C 65, 051305(R) (2002). 
Y. Yu and A. Bulgac, Phys. Rev. Lett. 90, 222501 (2003); 
ibid. Phys. Rev. Lett. 90, 161101 (2003); A. Bulgac and 
Y. Yu, Phys. Rev. Lett. 91, 190404 (2003); ibid. J. Low. 
Temp. Phys. 138, 741 (2005). 

The Many-Body X challenge problem formulated by G.F. 
Bertsch, see R.A. Bishop, Int. J. Mod. Phys. B 15, Hi, 
(2001). 

J. Carlson et al, Phys. Rev. Lett. 91, 050401 (2003). 
J. Carlson and S. Reddy, Phys. Rev. Lett. 95, 060401 

(2005) . 

O. Juillet, New J. Phys. 9, 163 (2007). 

J. Carlson et al, Phys. Rev. C 68, 025802 (2003); C. 

Lobo et al, Phys. Rev. Lett. 97, 200403 (2006). 

S.-Y. Chang and G.F. Bertsch, Phys. Rev. A 76, 

021603(R) (2007). 

J. von Stecher et al, arXi v:0705.0671V l. 

A. Bulgac and M. M. Forbes, Phys. Rev. A 75, 031605(R) 

(2007). 

R.G. Littlejohn and M. Cargo, J. Chem. Phys. 117, 27 



5 



(2002). 
[18] see online appendix. 

[19] A. Bhattacharyya and R.J. Furnstahl, Nucl. Phys. A747, 
268 (2005); ibid. Phys. Lett. B607, 259 (2005). 

[20] A. Bulgac et al, Phys. Rev. B 52, 16476 (1995). 

[21] Y. Nishida and D.T. Son, Phys. Rev. Lett. 97, 050403 
(2006). 

[22] J.E. Thomas et al. Phys. Rev. Lett. 95, 120402 (2005). 



Online appendix 



The emerging SLDA equations Eqs. ([8]) have been dis- 
cretized using the Bessel Discrete Variable Representa- 
tion (DVR) method outlined in Ref. 



171 ] . I have used 



two DVR bases, one with angular momentum I = for 
all even angular momenta and with I = 1 for odd an- 
gular momenta. In this representation only a part of 
the kinetic energy is a non-diagonal matrix, while both 
the meanfield and the pairing fields are diagonal. For 
angular momenta I > 1 either 1(1 + l)/2r 2 for even 
l's and [1(1 + 1) — 2]/2r 2 for odd I's was included into 
the potential. The number of discrete points was up 
to N max — 100, which is equivalent to the use of an 
harmonic oscillator basis with at least an equal num- 
ber of oscillator shells. Larger basis sets have been used 
as well. Comparable results were obtained with a ba- 
sis set of 35 points. However, unlike the decomposition 
of the quasi-particle wave functions in a harmonic oscil- 
lator basis, which has significant convergence problems, 
the use of Bessel DVR method is noticeable faster con- 
verging and for many quantities one can relatively eas- 
ily reach relative machine precision (sa l.Oe — 15). The 
iterative process was also significantly improved by the 
use of the Broyden method in determining the mean- 
field, pairing field and chemical potentials. A fully con- 
verged solution, at the level of 1.0e-9, was obtained in 
less than 20-30 iterations, irrespective of particle num- 
ber. The accuracy is typically improved by an order of 
magnitude every 2-3 iterations, once the process starts 
converging. All calculations were performed on a laptop, 
using a MATLAB program of about 400 lines, including 
the plotting and the comparison with GFMC/FN-DMC 
calculations. The part which solves the SLDA equations 
alone is about 200 lines. For a basis set of N max = 35 
points, a fully self-consistent solution for a given particle 
number is obtained in about 6 seconds, while for a basis 
set of N max = 75 points the time increases to about 30 
seconds and 55 seconds on average for N max — 100. Odd 
systems usually require more iterations. The program 
spends most of the time (more than 90 %) diagonaliz- 
ing the SLDA equations ([5]). A number of details can 
be found at www.phys.washington.edu/users/bulgac/ 
M edia/MM FAB.pdf. 

In order to describe systems with an odd particle num- 
ber one has to place a quasi-particle in a specific quan- 
tum state no and change the densities according to the 
following prescription: 



,(r)| S 



,(r)| 5 



(21) 



n(r)=2^K(r)| 2 + 

n 

r c (r) = 2]T |V V „(r)| 2 + |V M „ (r)| 2 - | Vv no (r)| 2 ,(22) 



"c(r) = 2j<(i>„(r) - u no (r)u* (r). 



(23) 
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It can be shown that this corresponds to exactly one par- 
ticle difference between the spin-up and spin-down parti- 
cle numbers. It also amounts to the fact that the single- 
particle quantum state n is excluded from the pairing 
mechanism. This is often referred to as blocking and it 
is the reason for the staggering in the ground state en- 
ergies of the even and odd systems. Essentially identical 
results can be obtained by introducing different chemical 
potentials for the spin-up and spin-down components. 

TABLE I: Table I. The energies E(N) calculated within 
the GFMC [lj], FN-DMC [H and SLDA. When two num- 
bers are present the first was calculated as the expectation 
value of the Hamiltonian/functional, while the second is the 
value obtained using the virial theorem, namely E(N) — 
mw 2 f d 3 rn(r)r 2 [II]. 
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The SLDA energies and the corresponding GFMC and 



FN-DMC results are presented in Table I for a DVR ba- 
sis set using N max = 100 points. In a calculation using 
N m ax = 35 points one observes differences typically at 
the level of 1.0e-4 or less for the energies, which can be 
ascribed to the weak dependence of the results on the 
cutoff energy E c . However, the agreement between the 
energy expectation values and the values calculated using 
the virial theorem degrades as N max is decreased. M.M. 
Forbes noticed that using a slightly smoothed cutoff pro- 
cedure eliminates a significant part of the energy cutoff 
dependence. For odd systems in particular we found that 
performing the calculations at a very small temperature 
T = 0.01 • ■ ■ O.lfhui eliminates the need to guess the opti- 
mal quantum numbers tiq of the extra quasi-particle. 
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FIG. 4: (Color online) The number densities for N = 8 (solid 
line), 14 (dashed line) and 20 (dot-dashed line), calculated 
within SLDA (blue line) and GFMC (red line), units as in Fig. 
[2] There are no shell effects in the SLDA number densities and 
their presence in the GFMC results is probably due partially 
to statistical reasons and/or partially to a less than optimal 
structure of the trial wave function. 

The number densities show a reasonable agreement in 
the surface region. However, the GFMC number densi- 
ties, unlike the SLDA number densities, show pronounced 
shell effects in the interior region. The magnitude of the 
pairing field for the smaller systems is such that it would 
not allow for an efficient pairing across harmonic oscilla- 
tor shells. The odd-even staggering of the total energies 
is however equally well manifest over the entire range of 
particle number. 



